A Multi-scale Approach for Simulations of Kelvin Probe Force Microscopy with 

Atomic Resolution 
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The distance dependence and atomic-scale contrast recently observed in nominal contact poten- 
tial difference (CPD) signals simultaneously recorded by KPFM using non-contact atomic force 
microscopy (NCAFM) on defect-free surfaces of insulating, as well as semiconducting samples, have 
stimulated theoretical attempts to explain such effects. Especially in the case of insulators, it is not 
quite clear how the applied bias voltage affects electrostatic forces acting on the atomic scale. We 
attack this problem in two steps. First, the electrostatics of the macroscopic tip-cantilever-sample 
system is treated by a finite-difference method on an adjustable nonuniform mesh. Then the re- 
sulting electric field under the tip apex is inserted into a series of atomistic wavelet-based density 
functional theory (DFT) calculations. Results are shown for a realistic neutral but reactive silicon 
nano-scale tip interacting with a NaCl(OOl) sample. Bias-dependent forces and resulting atomic 
displacements are computed to within an unprecedented accuracy. Theoretical expressions for am- 
plitude modulation (AM) and frequency modulation (FM) KPFM signals and for the corresponding 
local contact potential differences (LCPD) are obtained by combining the macroscopic and atomistic 
contributions to the electrostatic force component generated at the voltage modulation frequency, 
and evaluated for several tip oscillation amplitudes A up to 10 nm. For A = 0.1 A, the computed 
LCPD contrast is proportional to the slope of the atomistic force versus bias in the AM mode and to 
its derivative with respect to the tip-sample separation in the FM mode. Being essentially constant 
over a few Volts, this slope is the basic quantity which determines variations of the atomic-scale 
LCPD contrast. Already above A — 1 A, the LCPD contrasts in both modes exhibit almost the 
same spatial dependence as the slope. In the AM mode, this contrast is approximately proportional 
to A -1 ' 2 , but remains much weaker than the contrast in the FM mode, which drops somewhat faster 
as A is increased. These trends are a consequence of the macroscopic contributions to the KPFM 
signal, which are stronger in the AM-mode and especially important if the sample is an insulator 
even at sub-nanometer separations where atomic-scale contrast appears. 
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I. INTRODUCTION 



Kelvin probe force microscopy (KPFM), which was in- 
troduced twenty years ago,— £ has become an attractive 
non-contact technique to determine the electric surface 
characteristics of materials. Among many applications, 
this technique has been successfully applied for mapping 
local work function or surface potential variations along 
inhomogeneous surfaces of various materials^— For a 
conducting crystal, the work function corresponds to the 
energy difference between the vacuum level outside the 
surface at a distance large compared to the lattice spac- 
ing, yet small compared to the lateral dimensions of a 
homogeneous patch, and the bulk Fermi level. In this 
range, which is typical for conventional KPFM measure- 
ments, the potential acting on an electron approaches the 
local vacuum level and becomes constant, except in the 
vicinity of surface steps or patch boundaries. Differences 
between local vacuum levels are solely due to electro- 
static contributions which give rise to fringing electric 
fields around such boundaries. If the sample is covered 
by a thin overlayer of foreign material, the work func- 
tion can change owing to electron transfer and structural 
relaxation at the interface. 6 Similar changes can occur 
at the surface of a doped semiconductor, owing to band 
bending in a subsurface depletion layer. As long as elec- 



trochemical equilibrium occurs the Fermi level is aligned 
throughout the sample and with the Fermi level of the 
back-electrode. However, if the sample is a wide-bandgap 
insulator, e.g. an alkali halide crystal, this equilibration 
may require very long times, so that the bulk Fermi level 
is not well-defined. Charge rearrangements and relax- 
ation occur at the interface with the back electrode and 
cause an additive shift of the local vacuum level outside 
the surface with respect to the Fermi level of the back 
electrode. In a real, thick enough insulator with charged 
impurities, such a shift will also be affected by the dis- 
tribution of spatially separated charged defects at the 
interface, the surface and in the bulk of the sample, as 
well as by their slow diffusion over time^i 

When two separated conducting bodies, e.g. the probe 
tip of an Atomic Force Microscope (AFM) and the sam- 
ple, with different work functions <f>t and <p s are con- 
nected via back electrodes, electrons are transferred un- 
til the Fermi levels become aligned. The charged bod- 
ies then develop a contact potential difference (CPD), of 
Vcpd = {4>t—4>s)/e-, where e is the elementary charge and 
the sample is grounded. If the tip is biased at VJ, with re- 
spect to the sample, a finite electric field E ex V develops 
in the gap between them and causes an attractive elec- 
trostatic force proportional to V 2 where V — VJ, — Vcpd 
is their effective potential difference. If the sample is an 



insulator the same phenomenon occurs, but <j> a must be 
referred to the Fermi level of the back-electrode and is 
therefore affected by all the above-mentioned shifts, and 
so is Vcpd- It is then more appropriate to focus atten- 
tion on variations of Vcpd along the surface rather than 
on its absolute value which is affected by sample prepa- 
ration. 

In KPFM, a signal determined by this electrostatic 
force is compensated by applying a DC bias VJ, = Vcpd- 
For fast measurements the applied bias consists of an AC 
modulation voltage with angular frequency u> — 2nf in 
addition to the DC voltage: 



V b (t) = Vdc + Vac cos cot. 



(1) 



Assuming that the electric response is linear and in-phase 
with Vac , the electrostatic force acting on the tip can be 
decomposed into three spectral components: 



F(t) = F DC + F u cos cot + F 2u cos2wi. 



(2) 



The to component of the KPFM signal, which is propor- 
tional to (Vdc — Vcpd) Vac, is selectively detected by a 
lock- in amplifier and compensated by a feedback circuit. 
CPD variations along a surface can be conveniently 
measured together with its topography, 2 as determined 
by non- contact atomic force microscopy (NCAFM). In 
most state-of-the-art NCAFM experiments a micro- 
fabricated cantilever with a tip at its free end (typically 
etched out of doped single-crystal silicon) oscillates with 
a constant amplitude A at the frequency /i of a flexural 
resonance (usually the fundamental mode) 4^ Distance- 
dependent tip-sample forces cause a frequency shift A/i 
which can be very accurately measured using FM detec- 
tion (frequency demodulation)^ and used for distance 
control. In combined NCAFM-KPFM, the F^ component 
is simultaneously sensed; either the modulated deflection 
signal (Amplitude Modulation KPFM 11 ) or the modu- 
lation of the resonance frequency shift A/i (Frequency 
Modulation KPFMi 2 -) is actually detected.!^ In either 
case the amplitude of the signal at the modulation fre- 
quency / is proportional to (Vdc — Vcpd)Vac- Thus 
Vcpd can be recorded by continuously adjusting Vdc 
so that the modulation signal vanishes while scanning 
the tip parallel to the sample surface at a distance con- 
trolled by the (non-modulated) shift A/i. 9 Both modula- 
tion techniques are much faster and more sensitive than 
the direct method in which Vcpd is determined from 
the extremum of the parabolic A/i(T4) curve measured 
by slowly sweeping V\, at each measurement point liii^lS 
Potential artifacts of the modulation techniques^ are 
avoided in the direct quasistatic method. Because the 
scope of this article is primarily theoretical, we don't 
further consider such experimental difficulties, but fo- 
cus our attention on still controversial atomic-scale vari- 
ations of the so-called local CPD or Vlcpd on large 
defect-free surface areas. Thus we deliberately leave 
out local changes due to charged surface defects ^ 17 ' 20 
or adsorbate o 18 ' 21 which have recently attracted consid- 
erable attention, also in theor y 22 ' 23 



Atomic-scale variations of A/i can be detected by 
NCAFM on well-prepared surfaces in ultrahigh vacuum if 
the closest approach distance of the tip is smaller than the 
lattice spacing or the spacing between protruding atoms. 8 
The contrast in A/i then arises from short-range inter- 
atomic forces which begin to act in that distance range, 
while cantilever jump-to-contact is avoided if the total 
force remains much smaller than the maximum restoring 
force kA, k and A being respectively the flexural lever 
stiffness and oscillation amplitude. 9 Combined NCAFM- 
KPFM experiments have proven that FM-KPFM^^. 2 ! 
as well as AM-KPFM^ 8 .^ could detect lateral atomic- 
scale variations of Vlcpd in the range where A/i ex- 
hibits similar variations on surfaces of semiconductors, 
as well as of ionic crystals. Understanding the connec- 
tions between the observed contrast in Vlcpd and the 
atomic-scale variations of the electrostatic potential just 
outside the surface has been a challenging task, especially 
on unreconstructed cleavage faces of rocksalt-type crys- 
tals^ Above a flat homogeneous surface Vlcpd must, 
in principle, approach the corresponding Vcpd at some- 
what larger tip-sample separations. In practice, however, 
this ideal behavior is often masked by a slow dependence 
caused by the finite lateral resolution of surface inhomo- 
geneities, e.g. finite islands of materials with different 
work functions. This effect is less pronounced in FM- 
than in AM-KPFMJ^iii^ 2 . Several researchers devel- 
oped models and computational schemes based on clas- 
sical electrostatics which treated the tip and the sam- 
ple (sometimes also the cantilever) as macroscopic bod- 
ies in order to interpret the resolution of KPFM im- 
ages of inhomogeneous surfaces on lateral scales of sev- 
eral nanometers and abovei^ 3 . - — On the other hand, only 
few authors considered atomistic nano-scale tip-sample 
systems, either neglecting 1 -^ or including the macro- 
scopic contributions via simple approximations. In the 
first theoretical study of combined NCAFM-KPFM on 
an ionic crystal sample ^ 29 ' 42 a formally correct parti- 
tioning was proposed between capacitive and short-range 
electrostatic forces induced by the effective macroscopic 
bias V . This analytic treatment also provided qualitative 
insights into the origin of atomic-scale LCPD contrast, 
although underestimating the capacitive forces caused a 
quantitatively disagreement with experimental results as 
will be explained in subsection MI Al More reliable re- 
sults were obtained for a NaCl(OOI) sample interacting 
with a model tip consisting of a conducting sphere ter- 
minated by a NaCl cluster by allowing local atomic de- 
formations* 4 ^ These atomistic simulations were based on 
the SCIFI code 4 ^ which has provided detailed insights 
into NCAFM on ionic compounds i 45 ' 46 

In the present work, which is based on separate clas- 
sical electrostatics and ah initio calculations, we propose 
a more rigorous and accurate approach for coupling in- 
teractions acting on widely different length scales which 
leads to an unambiguous definition of Vlcpd- The out- 
line of this paper is as follows: in Section |TT] we discuss 
previous approaches, then present our own computation- 



ally simple, yet flexible finite-difference (FD) scheme with 
controlled accuracy to treat electrostatic tip-sample in- 
teractions on macro- and mesoscopic scales. Owing to 
electric field penetration into the dielectric sample, the 
tip shank and the cantilever significantly affect the ca- 
pacitive force and its gradient even at sub-nanometer 
tip-surface separations where atomic-scale contrast ap- 
pears. We also explain how the influence of the effective 
bias V can be included into atomistic calculations, as 
well as shortcomings of previous attempts to do so. In 
Section [IIII we critically discuss previous atomistic calcu- 
lations, as well as experimental evidence for short-range 
electrostatic interactions. Density functional calculations 
for nano-scale tip-sample systems are then discussed and 
illustrated for a realistic Si tip close to a NaCl(OOl) slab 
as an example of current interest. One important result 
is that the microscopic short-range force is proportional 
to V over a few volts; the corresponding slope is thus 
the basic quantity that should be extracted from KPFM 
measurements. In Section IIVI expressions for Vlcpd m 
AM- and FM-KPFM are obtained and evaluated, first 
for ultrasmall, then for finite tip oscillation amplitude 
A. Their magnitude and dependence on A are explained 
in detail in terms of the above-mentioned macroscopic 
contributions to the capacitive force. Experimental lim- 
itations and evidence for the predicted trends, as well 
as desirable measurements are also briefly discussed. Fi- 
nally, in Section [V] the main features of our approach 
and of our results are summarized, and conclusions are 
drawn. 



II. MACROSCOPIC ELECTROSTATIC 
INTERACTION 

A. Previous approaches 

Calculating the cantilever-tip-sample electrostatic in- 
teraction is, in fact, an intricate electrostatic boundary- 
value problem. The main difficulty is due to the pres- 
ence of several length scales determined by the nontrivial 
shape of AFM components, as well as to the distance- 
dependent redistribution of the surface charge density at 
constant bias voltage. In the case of conductive bodies 
with cylindrical symmetry, a simple assumption (uniform 
electric field along field lines approximated by circular 
arcs to their surfaces) led to an analytic expression for 
the force on the tip treated as cone with a spherical end 
cap4£ Recent numerical calculations^^ 8 - showed that 
Hudlet's expression is surprisingly accurate. Somewhat 
different analytical expressions and estimates for the lat- 
eral resolution in AM- and FM-KPFM were obtained for 
similar probes, also including a tilted cantilever^. For 
cylindrical geometries, many authors proposed numeri- 
cal schemes based on the image charge method which is 
applicable to simple geometries involving spherical and 
planar surfaces. 49 Thus Belaidi et al 50 placed N point 
charges on the symmetry axis and determined their po- 



sitions and strengths by forcing the potential on the tip 
surface to be V by a nonlinear least squares fit. The 
previously mentioned authors also described how con- 
tributions of the spherical cap, the tip shank and the 
cantilever to the macroscopic force lead to characteris- 
tic distance dependencies on scales determined by the 
geometry and dimensions of those parts. A linearized 
version of the numerical image charge method where the 
positions of axial point and line charges were fixed was 
applied to study tip-shape effects for conductive and di- 



electric sample. 



,37.51 



and thin films on conducting sub- 



strates^, also including the influence of the cantilever 53 . 
It is not known to what extent the boundary conditions 
must be satisfied for a given accuracy in the numerical 
image method, unlike in the analytic method where the 
positions and strengths of the image charges change with 
tip-sample separation and the boundary conditions are 
fully satisfied (see Appendix |BJ) . 

A more systematic approach to multi-length-scale 
problems is the boundary element method (BEM) 38 i 39 ' 48 . 
In this method the 3D (2D) differntial Poisson's equa- 
tion is transformed into 2D (ID) integral (Green's func- 
tions) equations on the surfaces of conductive or dielec- 
tric components, including CPD discontinuities and sur- 
face charges if desired 4°- The accuracy of BEM is con- 
trolled by the mesh resolution and is applicable to com- 
plex probe-sample systems, e.g. including a realistic 
cantilever 134 . The size of the resulting linear system of 
equations is small compared to volumetric discretization 
methods. However, because of the memory requirement 
of 0{N 2 ) to store the fully populated matrix and com- 
plexity of 0(N 3 ) to solve the linear equations, BEM has 
mostly been applied to systems with a relatively small 
number N of grid points, e.g. problems of high sym- 
metry and homogeneity for which it is feasible to de- 
rive the Green's function analytically. Somewhat earlier 
a few authors adapted Green's function methods devel- 
oped for more complex near-field optics problems to in- 
vestigate lateral resolution in KPFM on inhomogeneous 
samples^ 7 -. One advantage of BEM is that the LCPD 
of such samples can be expressed as a 2D convolution 
of the CPD and/or of a fixed surface charge distribu- 
tion with a point-spread function which depends only on 
the relative position of the scanning probei 38 ' 48 i 54 The 
distance-dependent lateral resolution can be quantified 
by the width of that function. Moreover, if one assumes 
that only one of those distribution is present, its can be 
determined by inversion of the BEM matrix upon dis- 
cretization on the adjustable BEM mesh. 40 

Conceptually more straightforward approaches involv- 
ing surface elements have been applied to conductive 
probe and sample systems. In the simplest one, the 
tip surface is approximated as a regular staircase (or, 
equivalently, as an array of capacitors in parallel) ) 13 ^ 4 ' 55 . 
More accurate methods rely on adjustable meshes. Thus 
the finite element method (FEM) was used to calculate 
the electrostatic force acting on a conical tip, 35 while a 
commercial FEM software was recently applied to sim- 



ulate a realistic cantilever and tip of actual shape and 
dimensions over a conducting flat sample with a CPD dis- 
continuity.— More sophisticated software packages have 
been used to solve the Poisson's equation in the presence 
of space charges, e.g. for structured samples involving 
doped semiconductor o 33 ! 57 . Numerical methods which 
involve 3D discretization require a very large number of 
grid points even if the mesh is carefully adjusted; the 
computational box must therefore be truncated at some 
finite extent. 



B. Implementation of finite-difference method 




As an alternative we present a finite-difference method 
(FDM) on a 3D non-uniform grid which is capable of 
dealing with realistic sizes of the cantilever, tip and sam- 
ple. Inhomogeneous metallic and dielectric samples as 
well as thin dielectric films on metal substrates, can be 
straightforwardly treated with this method. The most at- 
tractive feature of our FDM compared to FEM or BEM 
computations is its ease of implementation. Since the 
electrostatic potential varies smoothly and slowly at dis- 
tances far from the tip apex, we use a grid spacing which 
increases exponentially away from this region. Conse- 
quently, the number of grid points depends logarithmi- 
cally on the truncation lengths, and an extension of the 
computational box costs relatively few additional grid 
points. It allows us to simulate the cantilever as well as 
thick dielectric samples according to their actual sizes in 
experiments. 

The capacitance C(s) between the probe and the sam- 
ple back-electrode depends only on the tip-sample sepa- 
ration s, provided that their geometries are fixed/^ The 
macroscopic electrostatic energy due to the effective volt- 
age difference V = V& — Vcpd between the conducting 
tip and back-electrode is given by U c (s, V) — jC{s)V 2 . 
The electrostatic force exerted on the tip is proportional 
to the capacitance-gradient C'(s) 



W)~ff- 



§X=+^' 2 



(3) 



Similarly, the force-gradient is proportional to C"(s) = 
d 2 C/ds 2 . We wish to emphasize the difference between 
U and U c which leads to the positive sign on the RHS 
of Eq. (j3j) ; the reason is restated for convenience in Ap- 
pendix EJ The electrostatic energy 

U c {s,V) = ^J e(v)\V<S>\ 2 dv 

can be determined once the electrostatic potential 
<I>(r; s, V) is known at any point r in space. In general, 
when the dielectric constant e(r) varies in space, $ sat- 
isfies the generalized form of Poisson's equation 



FIG. 1. (color online) 2D (p, z) maps of the macroscopic elec- 
trostatic potential at three magnifications ( x 10 4 and x 10 5 in 
zoom-in windows). An effective bias V — 1 Volt is applied 
to the conducting probe while the back-electrode as well as 
the surrounding enclosure are grounded. The yellow region 
corresponding to $ = IV reflects the assumed cylindrically 
symmetric probe geometry: a cone with 15° half-angle termi- 
nated by a spherical cap of radius R — 20 nm. The cone is 
15nm high and attached to a disk of thickness 0.5 um. The 
radius of the disk is 35/xm which matches the area of a typical 
cantilever. The sample is a 1-mm thick dielectric slab with the 
relative permittivity e/eo = 5.9 of NaCl. The back electrode 
and the surrounding enclosure of height and radius 10 6 i? = 20 
mm are grounded (not shown). The sample- vacuum interface 
is indicated by the horizontal lines at z = 0. 



p being the charge density. Minimization of the energy- 
like functional 



I[tf(r)] 



f 



e(r) |V*| dv- / p^dr 



(5) 



subject to Dirichlct boundary conditions leads to <f>, the 
solution of the Poisson's equation Eq.(j4]) with the same 
boundary conditions^ Using a discretized variational 
approach, we therefore minimize the functional 



'({*=}) 






^n |Vtt£ - p„* n 



(6) 



V • [e(r)V$(r)] = -p(r) 



(4) 



On a non-uniform grid, we delimit the volume v n of the 
volume element assigned to node n by neighboring nodes. 
Then, ^ n , p n , e n and the electric field — V^n are eval- 
uated at the center of the volume element by linear in- 
terpolation between the nodes adjacent to n in orthog- 
onal directions. This ensures that the field is effectively 
evaluated to second order in the product of grid spac- 
ings and that discontinuities in V^n and e n at material 
interfaces are correctly treated. Although the formal- 
ism is general and can be applied to any 3D system on 
a judiciously chosen nonuniform 3D orthogonal grid, in 
the following examples we consider a cylindrically sym- 
metric setup without free charges in order to allow com- 
parison with most previous computations. In cylindrical 
coordinates, each volume element is a truncated tube of 



(z) 

height h k with inner and outer radii pi, Pi+i, respec- 
tively, and v n = 7r(p i+ i + pi)h\ p 'h!^\ h\ p > = p i+1 - pi 
Zk+i — Zk being respectively the radial and 



(*) 



and h 

vertical spacings of the appropriate nonuniform grid. The 
radial and vertical components of V\& are approximated 
on the circle of radius pi + O.bnf' at Zk + 0.5/i^ as 

(*<+!,», - *i,fc)//4 p) and (*i,*+i - *i,fc)/4 z) - Since the 

FD approximation of the electric field is a linear combina- 
tion of the potential values on nearest neighbor nodes, the 
functional in Eq. © is quadratic and the minimization 
condition dl /d^ n = yields a system of linear equa- 
tions A& = b where the vector b describes imposed 
boundary values and charge distributions. Because A 
is a sparse, symmetric and block-tridiagonal matrix, the 
system can be solved efficiently by an iterative procedure, 
which may, however, suffer from conditioning problems 
due to the nonuniformity of the grid. For an accurate 
solution, a mesh with high enough resolution is required 
in regions where $(r; s, V) varies strongly. We used the 
PARDISO package— to solve the resulting huge system of 
equations. An implementation of our FDM is distributed 
under GNU-GPL license as the CapSol code22. 

Once $(r, s,V=l) is determined for several separa- 
tions s, the system capacitance is obtained as C(s) = 
Je(r)|V<I>| 2 <ir ~ J2 n e " l v$ ln v a- Then a simple second 
order FD approximation is used to evaluate C'{s) and 
C"{s) from C(s). The electrostatic force acting on an 
arbitrary area S of a conducting part can also be evalu- 
ated as 



2eo 



o~(s) hdS, 



(7) 



where a(s) = —ed^/dn is the surface charge density 
guaranteeing that the tip surface is an equipotential, and 
n is the unit vector normal to the surface element dS. 
For a system with cylindrical symmetry the net force on 
a part of the probe delimited by two cylinders of radii 
px < P2 is vertical and given by F — 7re J P2 \V<&\ 2 pdp, 
however we prefer to use Eq. ([3|) to calculate the total 
macrosocopic force on the probe. In the following sub- 
sections we validate the performance of our FDM by com- 
parisons with previous results obtained by other methods 
for cylindrically symmetric systems. We mainly consider 
the macroscopic model system described in the caption 
of Fig. Q] which shows 2D (p, z) maps of the electrostatic 
potential computed by our FDM at three magnifications 
differing by five orders. The conducting probe consists of 
a conical tip terminated by a spherical cap of radius R 
attached to a cantilever modelled as a disk of the same 
area as a typical cantilever^ and the sample by a thick 
dielectric slab. Dirichlet boundary conditions are applied 
on a very large cylindrical box. Note that the grid spac- 
ing changes by six orders of magnitude (hundredths of nm 
around the tip apex to tens of pm near the box walls). 
The indented contours in the bottom right inset reveal 
the resolution of the finest grid, i.e. 0.02 nm. The con- 
tours in both magnified insets clearly show that for a 



separation of 1 nm a large fraction of the voltage drop 
occurs within the dielectric sample. Whereas the contour 
spacing between the tip apex and the surface is constant 
to a good approximation, it gradually increases inside the 
dielectric, in contrast to what occurs in a parallel plate 
capacitor. Actually the capacitance remains finite for an 
infinitely thick sample even in the (macroscopic) contact 
limit s — > (see Appendix IB"]) . 



C. Convergence and Accuracy 

Grid spacing: We first test our implementation for the 
problem of a conducting sphere of radius R separated by 
s from a semi-infinite dielectric surface for which an an- 
alytic solution of controllable accuracy is available (see 
Appendix [B| . A convergence analysis yields the param- 
eters needed to achieve a desired accuracy. Compared 
to the analytic solution, the convergence with respect to 
the finest grid spacing ho of the sphere-dielectric system 
(Fig. [2J shows a nearly quadratic error scaling for small 
separations as expected for a second order FDM. The 
accuracy could be improved by using a higher order ap- 
proximation for the electric field over further neighbor- 
ing points. However, a corresponding improvement of 
the approximation of curved surfaces on the orthogonal 
FD-mesh is then also required. Note that, for consis- 
tency, the surface of the sphere must be approximated as 
a staircase with variable step heights and widths which 
also change when the grid-spacing is changed. Thus the 
error scaling deviates somewhat from the ideal straight 
line and is no longer quadratic when the separation in- 
creases. The capacitance, force and force-gradient of this 
test system at a rather small separation of s = R/20 can 
be calculated within a relative error of 0.005 compared 
to the analytic solution if ho = i?/100. For larger sepa- 
rations, this accuracy is achieved even with a larger ho- 

For the cantilever-tip-sample system shown in Fig. Q] 
and the results in the next subsection, a uniform grid 
with h {p) = h,W = h = R/100 is used around the tip 
apex up to a distance of twice the tip apex radius in 
both radial and vertical directions. Outside this range 
the grid becomes gradually coarser with a growth factor 
of 1.01. In order to consistently preserve the shape of the 
tip approximated by the orthogonal mesh, the tip-sample 
separation must be changed in steps of ho- 

Space truncation: A convergence analysis with respect 
to the size of the computational cylinder is shown in 
Fig. [3] for the model system described in Fig. [TJ The 
capacitance approaches the same asymptotic value when 
the truncation length in a particular direction is increased 
while the other one is kept fixed and sufficiently large. If 
the computational box extends to 10 6 i? in the radial and 
vertical directions, the relative deviation of the capaci- 
tance from its asymptotic value is only 10 -7 (as indicated 
by the arrow in Fig. [3]) . We use these cutoff parameters 
in all subsequent FDM computations reported here. 
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FIG. 2. (color online) Convergence analysis with respect to 
the finest grid spacing ho for a conducting sphere of radius 
71 in front of a 1 mm thick dielectric of relative permittivity 
e/eo = 5.9. Points computed by our FDM for the macroscopic 
capacitance C, the force oc C" and force gradient oc C" are 
compared to the analytic solution for a semi-infinite dielectric 
described in Appendix [B] The sphere-surface separation is 
R/20 and the computational box extends to W e 'R in the radial 
and vertical directions. The straight line in the log-log plot 
indicates the expected quadratic error scaling. 
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FIG. 3. (color online) Convergence analysis with respect 
to the radial and vertical extents of the FDM computational 
box for the macroscopic system described in the caption of 
Fig. [U the tip-sample separation and finest mesh size being 
s — R/20 and ho = -R/100, respectively. The normalized 
capacitance of the system approaches the same asymptotic 
value upon increasing the truncation length in one direction 
while the other one is sufficiently large and fixed. Relative 
deviations with respect to the asymptotic value are shown in 
the inset. The arrow indicates the truncation length adopted 
in subsequent FDM computations. 



Comparison: In Fig. [4] we compare results obtained by 
our FDM with previous accurate BEM computations — 
for a system like in Fig. Q] but without the cantilever for 
a conducting and a dielectric (e/eo = 40) sample. The 
force and the force-gradient evaluated by the two meth- 
ods are in very good agreement for both kinds of samples. 
For the conducting sample, Hudlet's analytic approxima- 
tion^ deviates by only a few percent from the numerical 




FIG. 4. (color online) Normalized macroscopic electrostatic 
force (inset) and force-gradient computed by our FDM versus 
the normalized tip separation s/R from a dielectric (e/eo = 
40.0) and a conducting (e/eo = oo) sample compared to BEM 
computations (Ref.— ), as well as to Hudlet's approximation 
(Ref4£) in the second case (see text). The cantilever is ab- 
sent, as assumed in those two treatments, but the remaining 
parameters are as described in the caption of Fig. [T] 



results. In the following subsection we show that the con- 
tribution of the cantilever can be quite appreciable for a 
dielectric sample. 



D. Results 

The macroscopic electrostatic force and force-gradient 
versus the normalized tip-surface separation s/R for the 
system in Fig. Q] are shown in Fig. [5] for three differ- 
ent geometries: without, with a small and a large can- 
tilever modelled as disks of thickness 0.5 /im. The small 
disk radius is equal to the width of a typical rectangular 
AFM cantilever (20/im) while the total area of the large 
disk (of radius 35 /zm) matches the area of the rectan- 
gular cantilever. The presence of the cantilever increases 
the capacitance and the electrostatic force. Because the 
cantilever is more than 10 /im away from the surface, 
its contribution to the force is often considered constant 
for tip-sample separations smaller than R, and therefore 
does not contribute to the force gradient. Our calcu- 
lations [Fig. [5ja)] confirm that this is in fact true for a 
conductive sample. In this case, the main contribution to 
the force-gradient comes from the spherical cap, as can be 
seen from the solid line which corresponds to the analytic 
solution for a conducting spherical tip (see Appendix IB]) . 
However, the conical shank of the tip and the cantilever 
affect the force at large separations, as shown in the inset 
and noticed earlier i 33 i 34 ' 47 ' 50 On the other hand, if s/R is 
small, as shown in Fig. [5][b) and also emphasized in pre- 
vious work£i~— , over a thick dielectric sample both the 
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FIG. 5. (color online) Effect of the cantilever (size) on the 
macroscopic electrostatic force (inset) and force-gradient at 
different normalized tip separations from a conducting (a) and 
dielectric (b) sample. The cantilever is modelled as either a 
small or a large disk with radii of 20 and 35 /im, respectively. 
Other parameters are as in caption of Fig. [l] The solid lines 
show corresponding results for a tip approximated by a con- 
ducting sphere with radius R = 20 nm obtained by summing 
the analytic series for semi-infinite samples of both kinds (see 
Appendix |B)) . 



component of the electric field is two orders of magnitude 
stronger than the radial component parallel to the sur- 
face. These features are also clearly illustrated by the 
essentially equispaced horizontal equipotential contour 
lines in the vacuum region shown in the bottom inset of 
Fig. [TJ This important observation greatly simplifies the 
desired coupling to atomistic calculations: we can con- 
sider the electric field E z at the midpoint of the macro- 
scopic tip- surface distance s — z+h as a uniform external 
field acting on the isolated microscopic tip-sample system. 
The connection between those two scales is schematically 
illustrated in Fig. [7] 

Figure [5] shows that for a conducting sample the force 
gradient can be accurately described by a spherical tip if 
s < R, although the force itself is increasingly underesti- 
mated at larger separation a 15 ! 47 . In contrast, for a thick 
dielectric sample, the same description only provides the 
order of magnitude of Fm at small s/R, but exhibits a 
faster decrease with increasing separation and overesti- 
mates F' M . Figure IH] reveals that a spherical model tip 
overestimates the electric field E z under the tip at all 
separations, which then approaches V/R on the sphere 
(and zero on the surafce) when s 3> R. This occurs be- 
cause the induced surface charges can spread to the con- 
ical shank and the cantilever in the more realistic model. 
The contributions of those parts to the force Fm become 
nevertheless stronger than that of the sphere alone al- 
ready at small s/R. In general, if s/R — > 0, the electric 
field under the tip, hence the force and the force gradient 
are enhanced owing to an increasingly localized surface 
polarization of both tip and sample, but remain finite if 
the sample is a dielectric, as explicitly demonstrated by 
the solution for a spherical tip (see Appendix [Bjl . Com- 
parison with that solution (the solid curves in Fig. [5]) 
shows that even at small separations contributions from 
both the conical shank and the cantilever contribute to the 
force, whereas mainly the conical shank affects the force 
gradient. Hence, ignoring those contributions causes an 
overestimation of the force-gradient if the sample is an 
insulator. 



III. SHORT-RANGE ELECTROSTATIC 
FORCES 



force and the force- gradient are significantly decreased, 
owing to field penetration into the sample. 

A quantity of particular relevance in our multi-scale 
approach is the macroscopic electric field in the vacuum 
gap between the spherical tip end and the sample surface 
which polarizes the microscopic system. The variation of 
the electric field normalized to V/R at two points on 
the symmetry axis just below the tip and just above the 
surface is shown in Fig. [5] versus their normalized separa- 
tion. The same quantities shown magnified in the inset 
for nanotip separations relevant for atomic-scale contrast, 
i.e. z = s — h < 0.6 nm, differ little and drop only weakly 
with increasing z. In the same distance range the z- 



A. Evidence and previous models 

When an AFM tip approaches a surface, short-range 
forces contribute to the tip-sample interaction and give 
rise to atomic-scale contrast in NCAFM. Hereafter, F^ 
denotes the short-range force component perpendicular 
to the surface which can be extracted from measurements 
of A/i vs. the closest tip approach distance d in an os- 
cillation cyclej 61 i 62 The contrast observed in Vlcpd in 
the same distance range cannot only be due to the long- 
range electrostatic force, but must be due to a short- 
range bias-dependent force. Arai and Tomitori were the 
first to infer the existence of such a force from A/i(V&) 
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FIG. 6. (color online) Normalized macroscopic electric field 
just below the tip-apex and just above the dielectric sample 
surface (e/eo=5.9) versus their normalized separation for the 
probe described in the caption of Fig. [l] (curves with symbols) 
and for a tip approximated by a conducting sphere of the same 
radius (R = 20 nm) (continuous curves). Inset: zoom into the 
range where atomic-scale contrast appears; the electric field 
between the tip and the surface changes by only a few percent 
and is hence nearly uniform. 
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FIG. 7. (color online) Sketch of the AFM setup showing its 
macroscopic and microscopic parts on relevant scales. The 
macroscopic tip-sample separation is s — Z + h, where h is 
the nanotip height and z is the nominal distance (without re- 
laxation) between the apex atom and the surface. The macro- 
scopic electric field E obtained by solving Poisson's equation is 
applied as an external field to the atomistic subsystem shown 
in the zoom window. In the range z < 0.6 nm where atomic- 
scale contrast appears, the electric field in the vacuum gap can 
be considered to be uniform and equal to E z at the mid-point 
s/2 on the symmetry axis. 



curves recorded with a cleaned and sharpened silicon tip 
closer than 0.5 nm to a 7x7 reconstructed Si(lll) sur- 
faced In particular, above a Si adatom, they found a 
narrow peak growing with decreasing d superposed on 
the usual parabolic dependence around the plotted min- 
imum of — A/i(V(,) in their Fig. 1, i.e. for Vt, — Vcpd- 
Later the same authors pointed out that an even sharper 
peak appeared at the same bias in the simultaneously 



recorded tunneling current^ This seemingly supported 
their original suggestion that the additional attractive 
force causing the peak in — A/i(V;,) arose from the in- 
creased overlap due to the bias-induced energetic align- 
ment of dangling bonds states localized at the tip apex 
and on Si surface adatoms. The formation of a cova- 
lent bond between those states has been shown to be 
responsible for the observed NCAFM contrast on the 
7x7 reconstructed Si(lll) surfaced In extensive recent 
measurements on the same system, however, Sadewasser 
ct al reported parabolic Afi(Vb) curves, but detected a 
rapid drop by about - IV followed by a gradual increase 
in Vlcpd above a Si adatom with decreasing d in the 
range where the extracted short-range force showed a 
similar behavior J^ The apparent discrepancy with re- 
spect to Arai and Tomitori's observations is not so sur- 
prising because tunneling is seldom observed with clean 
silicon tips, although it is routinely measured in STM, as 
well as in NCAFM on conducting and even semiconduct- 
ing samples when using metal-coated silicon tipsj 66 i 67 
An appreciable position- and distance-dependent tunnel- 
ing current is, however, undesirable in dynamic atomic- 
scale LCPD measurements because it violates the basic 
assumption of an in-phase response to the AC voltage 
modulation by causing phase shifts which are difficult 
to compensate. This problem does not arise with insu- 
lating samples. Nevertheless, Arai and Tomitori's basic 
idea that bias-induced alignment of spatially localized 
surface states can lead to an enhanced site-dependent 
attractive force remains plausible even if a DC tunneling 
current cannot be sustained. Thus Krok and cowork- 
ers^ suggested that the lower LCPD which they found 
across protruding In rows on the c(2x8) reconstructed 
InSb(OOl) surface was due to a bias-induced local elec- 
tron transfer from a polar dangling bond on the elec- 
tronegative Sb atom presumably picked by the Si tip to 
the nearest electropositive surface In atoms. The same 
authors also showed that the LCPD contrast between 
different lateral positions decays exponentially with in- 
creasing d < lnm. 

Whereas bias-induced electron transfer is plausible 
for narrow-bandgap semiconductors like those previously 
mentioned, it is unlikely for overall neutral cleaved (001) 
surfaces of wide-bandgap insulators like alkali halides 
which neither have gap states, nor are reconstructed, but 
are only weakly rumpled. 68 In Rcf. 29 the atomic-scale 
LCPD contrast observed on KBr(OOl) was attributed to 
opposite surface cation and anion displacements in re- 
sponse to local electric fields induced by the macroscopic 
(in accordance with our definition) field. However, the 
authors approximated E z by the electric field V/R at 
the surface of an isolated conducting spherical tip, the 
local unit cell polarizability by the bulk crystal (Clausius- 
Mosotti) expression, and neglected the macroscopic sur- 
face polarization. Although essentially constant on the 
scale of a nanometer-size nanotip, the latter, together 
with E z is actually nonuniform on a lateral scale of order 
y/Rs for separations s<<-R. They evaluated the macro- 



scopic and microscopic surface charges densities a m and 
a^ induced on a conducting model tip by their E z and 
by the displaced surface ions, respectively. Using Eq. (|7J) 
they computed the modulation of the electrostatic force. 
After further justified approximations, they obtained op- 
posite LCPDs above cations and anions which increased 
exponentially with d. In a subsequent article^ 2 - the same 
authors added a macroscopic force roughly representing 
the interaction of the cantilever with the back electrode, 
but still obtained a surprisingly large maximum in the 
absolute LCPD for d ~ 0.6 nm. In a subsequent pub- 
lication^, even more reliable results were obtained for 
a cubic NaCl cluster partly embedded into a conduct- 
ing sphere interacted with a NaCl(OOl) sample similar to 
ours via empirical shell-model potentials. Cluster ions 
inside the sphere were fixed and allowed to interact with 
the protruding cluster ions which thus formed a nanotip 
with a net charge +e at the apex. The justification for 
such a model is that Si tips often pick up sample material 
and that simulations based on the same code produced 
reasonable results when compared to NCAFM measure- 
ments on alkali halides, e.g. on KBr(001)i 45 i 46 The re- 
sults obtained can be considered representative of what 
is expected for a strongly polar tip interacting with an 
ionic crystal. Recently a simpler model for such a tip 
(conducting sphere terminated by a point charge or a 
point dipole) could account for the observed variation of 
the LCPD over a few nanometers^ 

Earlier studies mentioned that the short-range tip- 
sample interaction is bias-dependent but provided no 
recipe to investigate it theoretically. Moreover, they 
did not clarify how long-range and short-range bias- 
dependent forces are connected and the role of each in 
the observed KPFM signals. In the following sections 
we answer all of these questions and obtain and analyze 
in detail theoretical expressions for the site-dependent 
LCPD. Our approach is not limited to particular materi- 
als, but results are presented for the system described in 
the following Section which is representative of a neutral, 
but polarizable reactive clean Si tip interacting with an 
ionic crystal. 



B. Density functional computations 

As illustrated in Fig. [7] our microscopic system consists 
of a nanotip of height h protruding from the spherical end 
of the macroscopic tip and of a wider two-layer slab of 
sample atoms. Computations are performed within the 
local-density approximation to density functional theory 
(DFT) using norm-conserving HGH pseudopotentials^ 
and the BigDFT package. 70 Relying on a wavelet basis 
set with locally adjustable resolution, this package calcu- 
lates the self-consistent electron density, the total energy 
and its electrostatic component with selectable bound- 
ary conditions^ i.e. periodic in two directions and free 
in the third in our case. This allows us apply an exter- 
nal field perpendicular to the surface without artifacts 



which can arise from periodic images in the z direction 
when using plane- wave of mixed basis sets. As already 
explained, the voltage biased macroscopic system deter- 
mines the uniform electric field E z oc V — V\, — Vcpd 
applied into the microscopic part (see Fig. [Jj . This pro- 
vides the desired well-defined relationship between the 
bias-voltage and short-range forces which was lacking in 
previous approaches to LCPD contrast based on DFT 
computations ! 16 ' 41 

Figure [5] illustrates the microscopic system used in the 
DFT computations reported here. The nanotip at the 
very end of a silicon tip is modelled as a cluster with 
a fixed (001) base of eight Si atoms with all dangling 
bonds passivated by H atoms in order to mimic the con- 
nection to the rest of the tip. The remaining Si atoms 
were prc-rclaxcd using the Minima Hopping Method 72 
previously employed to generate low-energy structures of 
silicon clusters and of similar model tips^ As in that 
work, the free Si atoms adopted a disordered configura- 
tion with several exposed under-coordinated atoms. In 
particular the protruding apex atom is threefold coor- 
dinated and hence has a dangling bond with a small 
dipole moment pointing towards the surface. As we ver- 
ified, a distance five times the lattice constant of NaCl 
is large enough to get rid of the electrostatic interaction 
between this nano-tip and its images in the periodic di- 
rections. Therefore our sample consists of two 10x10 
NaCl(OOl) layers containing 200 ions in total. For such a 
large system, we evidently perform calculations only at 
one single k-point, namely center of the surface Brillouin 
zone. With periodic boundary conditions applied along 
the main in-plane symmetry directions pre-relaxation of 
the sample only caused a small rumpling which preserved 
the basic periodicity of the truncated (001) surface. Al- 
though the silicon model tip and the sample were initially 
individually pre-relaxed, all tip and sample atoms were 
subsequently frozen in some of our KPFM simulations. 
In this way we could assess pure electronic polarization 
effects without effects due to the interaction-induced dis- 
placements of nuclei. 

The silicon model tip was positioned so that its fore- 
most atom was 0.65 nm above a sodium and chlorine 
surface ion, then moved towards the sample in steps of 
0.02 nm. At each step the Kohn-Sham equations are 
solved iteratively. The topmost layer of the Si tip to- 
gether with the passivating H atoms, as well as the bot- 
tom layer of the slab are kept fixed while other ions are 
free to relax until the Hellman-Feynman force exerted on 
each ion is less than 1 pN. This extremely tight toler- 
ance is required only when the relative variation of the 
force when the bias changes is very small. The force F^ 
exerted on the model tip is obtained by summing the 
z-components of the forces over atoms of the tip. Since 
the free atoms are well relaxed, their contribution to that 
force is not significant and was used as a measure of the 
error in forces. Figure [H] shows the microscopic force 
versus the tip-apex separation from CI and Na surface 
sites without applied electric field. The same procedure 
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FIG. 8. (color online) The microscopic Si-tip NaCl-slab sys- 
tem used in our ab initio DFT calculations. The apex of a 
silicon AFM tip is modelled as a pre-relaxed Si^gHts cluster. 
All eight atoms in the top (001) layer are passivated by hy- 
drogen atoms and kept fixed. The position of the foremost 
Si atom is (x,y, z), z being the height from the surface. The 
model sample consists of two NaCl(OOl) layers each contain- 
ing 10x10 ions with the bottom layer kept frozen. Periodic 
boundary conditions are applied along the x and y directions. 
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FIG. 9. (color online) Microscopic force on the Si nanotip 
above Na and CI surface ions from ab initio calculations with- 
out an applied electric field. Insets: variation of the force as 
a function of the macroscopic bias voltage at a tip-surface 
separation of 0.30 nm. The error-bar (not shown) is ~ lpN. 



is repeated at each tip-sample separation for a few field 
strengths E z determined as explained in subsection III Dl 
for effective biases — 2 < V = Vb — Vcpd < 2 Volts 
applied to the macroscopic tip. For such biases and dis- 
tances where F M becomes site-dependent, a nearly uni- 
form macroscopic electric field of ~0.15 V/nm occurs in 
the vacuum gap, as discussed in Sec. UU and illustrated 
in the inset of Fig. [5] No instabilites caused by elec- 
tronic and/or atomic rearrangements appeared in that 
range of parameters. The variation of the force at the 
particular separation z — 0.3 nm is shown in the insets 
in Fig. IHl In contrast to the macroscopic capacitive force, 
the short-range force depends linearly on the applied bias 
voltage. As explained elsewhere^ this linear term arises 
from the interaction between distance-dependent but V- 
independent net charge densities on the tip and sample 



surfaces with V~-induced changes on the opposite surface 
and with the macroscopic electric field. Earlier studies 
obtained such a term by treating native ions or charged 
atoms adsorbed on the sample surface and/or the tip 
apex as point chargesi 22 ' 23 ! 29 Deviations from the linear 
behavior could occur for larger biases, especially near 
instabilities, as observed in computations for a charged 
nanotip. 43 

The basic quantity which determines the deviation of 
the LCPD from the background CPD is the voltage- 
independent slope of the short-range force with respect 
to the applied voltage 

a(x, y, z) = gyF^x, y, z; E(V)). (8) 

As discussed in the Introduction, the background CPD 
is not a well-defined quantity for an insulator. For a real 
doped silicon tip-NaCl(OOl) sample, it would be different 
from the CPD of our microscopic system if charge equi- 
librium is achieved, as enforced by the self-consistency of 
the computations. Besides, no CPD is explicitly included 
in the description of the macroscopic system. Thus the 
effective bias V = Vb — Vcpd would differ from that in 
a real system. Nevertheless, as long as this bias is in the 
Volt range, the slope a is unaffected. 

Fig. [TUTa) shows that the slope a exhibits a character- 
istic site-dependent distance dependence at separations 
less than 0.5 nm, and is larger above the more polarizable 
CI ion. The underlying physics will be explained else- 
where^ The microscopic force-gradient F' is also a lin- 
ear function of bias voltage, and the distance-dependence 
of its slope a'(x,y,z) — dF'/dV, approximated to sec- 
ond order by linear interpolation between adjacent points 
on both sides of a given z- value, is shown in Fig. [TUtc). 
Figures [TUfb) and [TUH) show that a and a' are weaker 
if relaxation is allowed but that contrast appears below 
nearly the same distance and exhibits almost the same 
distance dependence. Thus, for the assumed neutral Si 
nanotip, the contrast is mainly due to electronic polar- 
ization rather than to bias-induced ion displacements. 

In the approximation that the macro- and microscopic 
systems are coupled only through the macroscopic elec- 
tric field, the z-component of the total force exerted on 
the tip is 



F = F M (s; V) + F^x, y, z; V) + F vdW (s) 



(9) 



where s = z + h and V = Vb — Vcpd- The long-range 
van der Waals force F vc tw is bias- and site- independent, 
being only a function of the mesoscopic geometry and 
is therefore henceforth ignored, although it affects the 
resonance frequency shift A/ in a NCAFM measurement. 
The macroscopic force Fm is capacitive (oc V 2 ) while the 
microscopic force F^ has been shown to be linear in V . 

Three additional corrections couple the bias-dependent 
macro- and microscopic forces. The first correction 
5CV 2 /2 is due to an additional capacitive contribution 
caused by the presence of a polarizable nanoscale object 
in the gap between the macroscopic bodies. Owing to 
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FIG. 10. (color online) Distance-dependence of the slopes 
a = dF^/dV and a = dF'^/dV above Na and CI surface ions 
with (a,c) and without(b,d) relaxation of the free atoms and 
ions during tip approach. The difference (contrast) between 
Na and CI sites is shown by red (filled) symbols. 



the small lateral dimensions of the nanotip compared to 
the radius of the macroscopic tip end, this correction is 
small, 23,74 although it can become appreciable and site- 
dependent if the nanotip apex is charged. 43 

The second correction arises only in that case or if the 
nanotip has a large net dipole moment 75 . This leads to a 
site-independent LCPD with an approximate power-law 
approach towards a background CPD of several Volts. 
The interaction of the nanotip charge distribution with 
the macroscopic field E could in principle be included in 
our description at separations s where E can no longer be 
considered uniform. In that range, however, the charge 
or dipole might be approximated as point objects, as jus- 
tified in the case of a conducting sample in the Sup- 
plemetary Material of Ref. 22 . Because the charge or 
dipole are intrinsic, the interaction is proportional to V, 
so that this correction would give rise to long-range con- 
tributions to the slopes a and a'j 22 ' 23 In the case of our 
neutral Si nanotip, this correction is small. 

The third correction arises because in reality the nan- 
otip is in electrical contact with the macroscopic tip, so 
that the electron density at the interface differs from 
that near the top of our isolated silicon cluster. How- 
ever, this model tip is large enough, so that the charge 
distribution near the apex, which dominates F^ is not 
much affected, in contrast to models with smaller model 
tips. The self-consistently determined microscopic elec- 
tric field between the nanotip apex and the sample sur- 
face differs from the applied macroscopic field E z , but 
this effect is already included in the computed micro- 
scopic force. 
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FIG. 11. (color online) (a) Sketch of the cantilever oscillating 
in its fundamental mode of a tip with finite amplitude A at 
the end . The weight functions used for calculating the cycle 
averages in Eqs. (|16I17[) are shown as functions of £ = z — d— A 
where d is the closest tip apex-sample separation. Dependen- 
cies of the first (b) and second (c) spatial derivatives of the 
capacitance on the macroscopic separation s = z + h calcu- 
lated for the setup in Fig. [TJ and of their cycle averages (d,e) 
on d for tip oscillation amplitudes A =0.01, 0.1, 1 and 10 nm. 



IV. AM AND FM KPFM SIGNALS AND LOCAL 
CONTACT POTENTIAL DIFFERENCES 



A. Ultrasmall amplitude limit 

The force gradient is more sensitive than the force to 
short-range interactions which are responsible for atomic- 
scale contrast in NCAFM and KPFM. Direct detection 
of the gradient is in principle possible if the variation of 
F^ over the peak-to-peak oscillation amplitude is linear, 
e.g. if 2A is comparable to the spacing 0.02 nm of the 
computed points in Fig. We first consider this simple 
limit which is commonly assumed in the KPFM litera- 
ture, but is seldom achieved in NCAFM experiments. As 
explained in the Introduction, Vlcpd is operationally de- 
fined by nulling the KPFM signal generated by the force 
component F u at the modulation frequency. Assum- 
ing that the response Vac is linear and instantaneous, 
Fui = (dF/dVb)VAc, an d the deflection signal detected in 
AM-KPFM would be proportional to 



F„ 



[C'(z + h) {V DC - Vcpd) + a{x, y, z)} Vac (10) 
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FIG. 12. (color online) Calculated deviations AVlcpd for 
AM- (left column) and FM-KPFM (right column) versus clos- 
est tip apex-sample for tip oscillation amplitudes A = 0.01 
nm (a.,e),A — 0.1 nm (b,f), A — 1 nm (c,g), and A = 1.0 nm 
(d,h). In (e,f) the dashed horizontal lines indicate the range 
of validity of our DFT calculations (±2V) . 



longer hold in the case of a conductive sample or thin 
dielectric film on a conductive substrate. 

In FM-KPFM the contribution of the modulated elec- 
trostatic force component F u to the frequency shift of 
the first resonant mode A/i is detected and nulled. In 
the ultrasmall amplitude limit A/i is proportional to the 
force-gradient 1° and would therefore be nulled if 

K = [C"(z + h) (V DC - V CPD ) + a'(x, y, z)\ Vac = 0. 
The FM-counterpart of Eq. (fP2")l is therefore 



av fm 



LCPD 



a'(x,y,z) 
C"(z + h)' 



(13) 



Again, as illustrated by by Fig. [TUtc) and by the points 
for ^4=0.01 nm in Fig.fTTTe). the site- and distance depen- 
dence of this deviation is determined by a'(x, y, z), but 
AV£cp D again differs from the numerator by an almost 
z-independent factor. The calculated LCPD deviations 
for A = 0.01 nm in the AM and FM modes are plotted in 
Figs.Q^a) and (e). Note that AV£^p D would be about 
hundred times stronger and would exceed the range of 
validity (±2 V) of our DFT computations (limited be- 
tween the horizontal lines in Figs. [T2T e.f)). as well as the 
measured results in experiment, hence cannot be trusted. 
Therefore it is important to consider averaging over the 
range covered by the finite tip oscillation. 



B. Finite amplitude expressions 

In NCAFM with cantilevers the oscillation amplitude 
A is between several and a few tens of nanometers, so that 
the macroscopic capacitive electrostatic force can change 
by several orders of magnitude over an oscillation cycle. 
In practice, the detected AM and FM KPFM signals are 
given by differently weighted averages, namelj*2£ 



in the ultrasmall amplitude limit, and would be nulled if 

a(x,y,z) 



Vdc = Vcpd — 



C'{z + h)' 



(11) 



Because the background Vcpd is not well-defined, and 
only a(x, y, z) is site-dependent, we consider only the de- 
viation of Vlcpd from Vcpd which is responsible for 
atomic-scale contrast, i.e. 



AV AM 
n - v LCPD 



a{x,y,z) 
C'(z + hY 



(12) 



For a thick dielectric sample, as illustrated by Fig. [POT a) 
and by the points for ^4=0. 01 nm in Fig. mT d). the z- 
dependence of C" is weak over the range where a(x, y, z) 
is appreciable (s = z + h < 1 nm). Therefore AVlcpd 
differs from a[x, y, z) by an essentially z-independent fac- 
tor. Depending on the nanotip height h, this may no 



1 f 2 " 
(F u ) =— / F u [d+A(l + coa 



and 76 



kA 



/ 



2W 



F ul [d + A(l + cos ^)] cos 



where k is the flexural stiffness of the cantilever and 
d — z min is the closest tip apex-sample separation. Sub- 
stituting the force from Eq. flTtfl) and setting these aver- 
ages to zero, one obtains 



AV AM 

AAV LCPD 



AV FM 

LAV LCPD 



(g(x,y,z)) 

'{C'{z + h)Y 

(g'(x,y,z)) 

(c»(z + h)y 



(14) 
(15) 
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where the cycle averages defined as 

(g) = - f W(()g(d + A + OdC 



(16) 



(flOs-L / (W(Og(d + A + Od( 
kA J_ A 



irA 2 



y/A*=(*g'(d + A + QdC (17) 



depend both on d and A. As depicted in Fig. ITlTa). 
C = z-(A+d) whereas W(Q = (A 2 ~C 2 )~ 1/2 is a weight 
function with square root singularities at the turning 
points of the oscillation. The expression on the second 
line of Eq. (|17p justifies the notation (g 1 ) and shows that 
this quantity tends to g'(d+A) when A — > 0, besides help- 
ing to relate the distance dependence of AVlcpd to those 
of a'(x, y, z) and C"(z + h). However, because a(x, y, z) 
is computed with high precision, whereas a'(x, y, z) is ob- 
tained by interpolation, we use the expression on the first 
line for numerical purposes. Furthermore, since a(x, y, z) 
is known only at equispaced separations z% where the 
DFT computations have been performed, the integrals 
in Eqs. (|16H17[) must be discretized. The adopted pro- 
cedure, which deals with the singularities of the weight 
function W(£) at the integration limits^ is presented in 
Appendix [Cj There we also show that the discretized ver- 
sion of the expression in the first line of Eq. (fl7|) reduces 
to the second order FD approximation of g'(d + A) when 
2A matches the spacing between adjacent z% values, in 
accordance with the expression on the second line. 



C. Results 

Owing to the very different z-dependencies of a{z) and 
C'(z + h), shown respectively in Figs. ITUTa) andQjJb), 
their cycle averages depend in different ways on d and 
A. The same holds for a'(z) and C"[z + h), shown re- 
spectively in Figs. HUT c) and [TTT cV Figures flTT d) and 
\TT\ e) show the cycle averages of C and C" versus the 
closest tip-apex approach distance d for oscillation am- 
plitudes A = 0.01, 0.1, 1 and 10 nm, whereas the cycle- 
averages of Vlcpd calculated from Eqs. (|14I15|) are plot- 
ted in Fig. QJ] for AM-KPFM (left column) and FM- 
KPFM (right column) for the same amplitudes in the 
range where a(z) is finite. In that range, the cycle aver- 
ages for A = 0.01 nm agree with the non-averaged quan- 
tities. Since the primary quantities were calculated at 
points spaced by 0.02 nm, this is not surprising in view 
of the remarks at the end of the preceding subsection. 
Thus, apart from small deviations introduced by the dis- 
cretization procedure, the points in Figs. IT27 a) and lT27 e) 
which were actually calculated for A = 0.01 nm coincide 
with those given by Eqs. (|12ll3p , and exhibit essentially 
the same distance dependencies as a(d) and a'(d), as al- 
ready discussed in the subsection IIV Al 

Already above ^4=0.1 nm, however, the LCPD con- 
trasts in both modes exhibit almost the same spatial de- 



pendence as a(d), although their respective magnitudes 
decrease if A is increased. Nevertheless, AV^p D sig- 
nificantly exceeds AV^pp,; this can be understood as 
follows. As seen in Figs. [TTT d) and [TTT e). (C") drops 
much faster than -(C) if A is increased. As explained in 
the discussion of Fig. [HJb) this behavior reflects the in- 
creasing influence of the relative contributions of the tip 
shank and of the cantilever to C'(z + h) in the range cov- 
ered by the peak-to-peak oscillation. Especially (C) is 
affected by the cantilever contribution which causes the 
very gradual levelling of C'(z + h) apparent in Fig.QjJb). 
As seen in Fig. Illf c). this slowly varying contribution 
tends to cancel out in C"(z + h), and, according to the 
second line in Eq. (fTT|) , in (C") as well. 

On the other hand, (a) and A(a') essentially coincide 
once a exceeds the range where a is noticeable. Indeed, 
the main contributions to those averages come from the 
vicinity of z — d where the integrands in Eqs. (|16|) 
and (fT7|) (first line) match. Expanding W(() about 
this turning point, one finds that (a) ~ A~ x l 2 whereas 
(a') - A- 3 ' 2 , just like A/i behaves in NCAFM^ Ac- 
cording to Fig. Illf b.c) the same argument cannot be ap- 
plied to (C) for A < 10 nm, and not at all to (C) 
because C'(s) varies only slowly up to s = R = 20 nm. 
Fig. Q2] shows how the finite oscillation amplitude affects 
the relevant cycle averages, as well as AVlcpd in the AM 
mode (left column) and in the FM mode (right column) 
at the closest tip apex-sample separation d = 0.30 nm 
indicated by arrows in Fig. [5] 

The same trends persist at all separations d < 0.5 nm 
where LCPD contrast appears, (a) drops as A~ x / 2 , and 
(a 1 ) drops as A~ 3 / 2 already beyond A — 0.1 nm, while 
(C) varies only little and (C) begins to drop somewhat 
slower than A -1 only above A = 1 nm. The resulting 
amplitude dependencies in both modes reflect the differ- 
ent dependencies of the numerators and denominators in 
Eqs. (fT4lfT5l) . 



D. Discussion and Experimental Limitations 

Expressions formally similar to Eqs. (|14H15p were ob- 
tained by Nony et al^ who also noticed that (a) and 
A(a') almost coincide when A exceeds a few nanometers. 
However, their denominators came from a short-range 
polarization contribution ex V 2 to the microscopic force 
F^ rather than from the much larger capacitive force 
F M . This results in a comparable AVlcpd for AM and 
FM modes if A exceeds a few nanometers. However, by 
including the correct Fm and taking into account the 
different amplitude dependencies of the denominators in 
Eqs. p4j|15[) . we conclude that the contrast should re- 
main larger in the FM than in the AM mode for a given 
closest approach distance d and oscillation amplitude A. 
This prediction is independent of the particular system 
considered, but the mode-dependent signal to noise ratio 
must also be considered. Thus Kawai et al^ calculated 
the minimum detectable CPD as a function of A and 
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FIG. 13. (color online) Amplitude dependencies of the cycle 
averages (o) and (a') (a,b), (C') and (C") (c,d) and of the 
resulting deviations AV L cpd an d AV L cpd ( e if) at a closest 
tip apex separation of d=0.3 nm above CI and Na surface 
sites. 



showed that it is smaller in the AM mode. Taking into 
account the discussions of Figs. [5l and [Til (C) would be 
larger if the cantilever area is larger whereas (C") would 
be unaffected, whereas both quantities would be larger if 
the cone angle is broader or if the sample is a metal rather 
than an insulator, but (C") would be more strongly af- 
fected. On the other hand (a) and (a') would be larger 
if the tip apex is charged^ rather than neutral, or if the 
sample is a semiconductor with a reconstructed surface 
which exposes partially charged species like Si(lll) 7x7 
16 ' 30 . From this point of view the system studied here 
is especially challenging. Furthermore, the contrast ratio 
slowly decreases if A is increased, e.g. by a factor which 
drops from about 100 to 10 for oscillation amplitudes 
between 0.01 and 10 nm in our example. 

For a meaningful comparison with NCAFM-KPFM 
measurements it is important to take experimental lim- 
itations into account. In view of the long-range LCPD 
variations due to surface and bulk inhomogeneities on 
real samples, one should compare computed atomic-scale 
LCPD variations with the difference between the LCPD 
measured at sub-nanometer separations d in the middle 
of a flat homogeneous island or terrace and the extrap- 
olated long-range, essentially site-independent LCPD. 
This procedure would also suppress most of the long- 
range contributions to (a) and (a 1 ) which would arise 
in the case of a charged or strongly polar tip22. More- 



over, the comparison should be done with the same tip 
at constant d (slow distance control) because atomic- 
scale variations of d at constant Afi(x,y,d) would in- 
duce such variations in the LCPD even if the latter is 
site-independent but has a different distance dependence 
as A/i. 

For the distance controller to function properly, A/i 
must be chosen on the branch where this frequency shift 
becomes more negative if d is decreased. Furthermore, 
the maximum restoring force kA must be much larger 
than the maximum tip-sample attraction^. For mea- 
surements with standard NCAFM cantilevers (k ~ 20- 
40 N/m) this criterion is typically satisfied by using os- 
cillation amplitudes A > 5 nm, and atomically resolved 
imaging is typically performed at distances d ~ 0.4- 
0.5 nm. According to Fig. [T3] the LCPD contrast which 
is then predicted to be 20-100 mV in the FM mode and 
a few mV in the AM mode approaches the experimental 
limits in both modes. Even if the AM-KPFM signal is en- 
hanced by setting the modulation frequency at the second 
fiexural resonance of the cantilever, the LCPD contrast 
predicted by our model would remain the same. This 
contrast would be stronger if the tip were charged. Unfor- 
tunately, available data showing atomic-scale contrast on 
(001) surfaces of NaCl and KBr is insufficient for a mean- 
ingful comparison between AM and FM KPFM. How- 
ever, LCPD maps obtained with sputter-cleaned Si tips 
and similar measurement parameters on Si(lll) 7x 7 sur- 
faces show that the contrast between Si adatoms and cor- 
ner holes in the FM-modoi^ is about ten times stronger 
than in the AM-mode^. Moreover, data obtained from a 
direct determination of the maximum of A/i versus bias 
voltage Vb agreed well with those obtained by nulling the 
FM-KPFM signal at the modulation frequency^. 

The sizable LCPD contrast of several Volts predicted 
in the FM mode for amplitudes A < 0.1 nm should, how- 
ever, be readily observable when using a tuning fork in- 
stead of a cantilever. Owing to the much higher stiff- 
ness k ~ 1800 N/m of this deflection sensor, the above- 
mentioned criterion can be satisfied with such ampli- 
tudes close to the ultrasmall limit 9 . Combined NCAFM- 
KPFM measurements using such tuning forks with Ptlr 
tips have only been done at low temperature by the time- 
consuming direct method mentioned in the Introduc- 
tion! 17 ' 18 Unfortunately, no FM-KPFM measurements 
showing atomic-scale LCPD contrast on alkali halide 
(001) surfaces have so far been reported. 



V. SUMMARY AND OUTLOOK 

We proposed a general multiscale approach to compute 
electrostatic forces responsible for atomic-scale contrast 
in KPFM performed simultaneously with NCAFM. Al- 
though attention is focused on insulating samples and 
results are presented for a silicon tip interacting with 
a NaCl(OOl) sample, the approach is not restricted to 
particular sample or tip materials. The problem is split 
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into two parts coupled in a remarkably simple but novel 
fashion. First the electrostatic problem of the voltage- 
biased AFM probe (including the tip and the cantilever) 
against the grounded sample, treated as macroscopic 
perfect conductors or insulators, is solved by a finite- 
difference method with controlled accuracy on a non- 
uniform mesh. The method is capable of treating com- 
plex geometries with widely different dimensions, but is 
illustrated here for systems with cylindrical symmetry. 
The solution yields the electric potential and field dis- 
tributions and the capacitance C(s) of the system from 
which the electrostatic force Fm acting on the probe and 
its gradient are calculated as functions of the macroscopic 
tip-sample separation s. By comparing results obtained 
with and without the cantilever, as well as with the an- 
alytic solution for a tip approximated by a conducting 
sphere in Appendix[B] the contributions of the cantilever, 
the conical tip shank and of its spherical end can be rec- 
ognized. If the sample is a thick insulator, all three af- 
fect the macroscopic force, whereas the last two affect the 
force gradient even at sub-nanometer separations relevant 
for atomic-scale contrast. 

Instead of the bias voltage V>, the nearly uniform elec- 
tric field obtained in that range is then applied as an ex- 
ternal field to the microscopic part which can be treated 
by empirical atomistic or first principles methods. The 
ab initio BigDFT wavelet code employed here enables 
one to compute the short-range bias-dependent force on 
the tip apex represented by a cluster with an unprece- 
dented accuracy of 1 pN per atom. For the Si-nanotip- 
NaCl(OOl) system considered here, this microscopic force 
F/j, is a linear function of the bias in the investigated range 
Vb — Vcpd = ±2 Volts. We argue that this is a general 
result, except close to atomic-scale instabilities caused 
by strong enough forces which could arise at very small 
separations and/or very large effective biases. 

Adding the macroscopic and microscopic bias- 
dependent forces, expressions are obtained for the KPFM 
signals in the AM and the FM modes. The atomic-scale 
deviation AVlcpd of the local CPD from its common 
value at large separations is the ratio of the derivatives 
a = dFft/dVb and dC/ds averaged over the tip oscilla- 
tion amplitude with different weights in AM- and FM- 
KPFM, as described by Eqs. (I14JI15J) . On the other hand, 
we explain the amplitude dependence of the atomic- 
scale LCPD contrast in both modes and predict that 
for typical amplitudes used in measurements with stan- 
dard NCAFM cantilevers, this contrast should be much 
stronger in the FM mode. This is a consequence of the 
contributions of the cantilever and the tip shank to the 
KPFM signal in the AM mode, which are stronger on in- 
sulating samples. The same conclusion has previously 
been reached in comparisons of AM- and FM-KPFM 
measurements of long-range LCPD variations; such vari- 
ations are caused by interactions of the biased probe 
with CPD inhomogeneities and surface charges on scales 
of several nanometers and above on conducting samples 
partly covered with ultrathin overlayers of different ma- 



terials 13 ' 31 . However, the strong mode-dependent influ- 
ence of distant contributions to dC/ds on the atomic- 
scale LCPD contrast has, to our knowledge, not been 
recognized because previous work on this topic assumed 
that only the tip apex mattered at sub-nanometer sepa- 
rations. 

Because AVlcpd depends on measurement param- 
eters, it is desirable to extract the more fundamental 
quantity a — dF^/dVb from combined KPFM measure- 
ments, just like the microscopic force F^ is extracted 
from NCAFM measurements using, e.g. a widely ac- 
cepted inversion algorithm 62 or one based on the direct 
inversion of the discretized version of the first line of 
Eq. flTS)) described in Appendix [C] by back-substitution 77 . 
Since AVlcpd is predicted to be stronger in FM-KPFM, 
whereas its distance dependence is governed by the 
weighted average (a') modes, the most appealing way 
to obtain a(d) would be to extract a' then integrate it 
from the range where AVlcpd vanishes down to the de- 
sired separation d. The averages (a') and (C") can be 
separately obtained from direct measurements of the fre- 
quency shift A/i as a function of bias — , namely from 
the shift of the maximum and the curvature of parabolic 
fits at several (x, y, d) positions. The signal/noise ratio 
of those averages can be improved by using AC modula- 
tion and lock-in detection at the modulation frequency. 
The averages could then be determined from the zero in- 
tercept Vlcpd and the sl °P e of the FM-KPFM signal 
(A/ w ) versus DC bias. An analogous procedure could be 
applied to determine (a) and (C) from the AM-KPFM 
signal (F u ), then a itself by inversion, using suitably 
modified algorithm s 77 ' 78 . Because the AM-KPFM sig- 
nal/ratio is much superior if the modulation frequency / 
is at the second cantilever resonance 30 , AV/^rd could 
be determined more accurately even if it is smaller than 
in FM-KPFM. In any case, note that the slope a reflects 
variations of the electrostatic potential outside the sam- 
ple surface which are, however, locally enhanced by the 
proximity of the tip apex. Since the latter is in turn also 
polarized and deformed—, a cannot simply be described 
as the convolution of the unperturbed electrostatic po- 
tential with a merely distance-dependent tip point-spread 
function, as in macroscopic electrostatics^ 

Complications due to averaging over the tip oscilla- 
tion amplitude are to a certain extent avoided with tun- 
ing fork deflection sensors which enable direct measure- 
ments of (A/ w ) vs. bias, using amplitudes approaching 
the ultrasmall limi t 17 ' 18 . Spectacular results have thus 
been obtained on isolated molecules adsorbed on a thin 
epitaxial NaCl(OOl) film by using tips with well-defined 
apex species stable at low temperature 2 ^. Most recently, 
Vlcpd contrast reflecting changes in the intramolecular 
charge distribution has been observed upon a configura- 
tional switch triggered by a judiciously applied pulse^ -. 
Our results shown in Figs. [T27 e), [T2lf) and [T"3ld) show 
that AV£cpd ano - a ' s ^ m bave a significant amplitude 
dependence between A = 0.1 and 0.01 nm, so that inver- 
sion is still necessary to obtain accurate results for typical 
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amplitudes used with tuning fork sensors. 

Since such measurements use hard metal tips, while 
metal-coated tips are also used in NCAFM and/or 
KPFM measurements with cantilevers it would desirable 
to develop appropriate nanotip models and to perform 
simulations like those described here. In particular, the 
recently fabricated sharp and stable W and Cr coated 
silicon tip o 67 ' 75 and the stable atomic-scale resolution 
achieved with Cr-coated cantilevers at separations ex- 
ceeding the usual range d < 0.5 nm merit further at- 
tention. Intentionally picked atoms or molecules at the 
apex would be worth studying in a further step. An- 
other class of systems which merit further investigations 
involve silicon nanotips with a picked-up cluster of for- 
eign material, NaCl in particular, which have so far been 
studied by DFT in the absence of a sample 81 or repre- 
sented by a cluster of the same material as the sample 
using empirical interaction potential o 43 ^ 45 . 

Note finally that all macroscopic probe models, includ- 
ing ours, provide a better description of metallic or metal- 
coated tips than of real silicon tips. Indeed, even if the 
native oxide is removed by sputtering, a silicon layer of 
few nanometers depleted of charge carriers still separates 
the tip surface from the highly doped conducting tip in- 
terior. Although it was taken into account in previous 
treatments of KPFM of semiconductor devices^, this de- 
pletion layer remains to be included when modelling Si 
tips, e.g. by allowing a smaller effective radius R of the 
equipotential at the applied bias voltage and a larger ef- 
fective separation s from the sample surface. 



ACKNOWLEDGMENTS 

This work has been supported by the Swiss National 
Science Foundation (SNF) and the Swiss National Cen- 
ter of Competence in Research (NCCR) on Nanoscale 
Science. The CPU intensive computations were done at 
the Swiss National Supercomputing Center (CSCS) in 
Manno. 



Appendix A: Sign of the macroscopic electrostatic 
force 



Using the virtual work method, the macroscopic elec- 
trostatic tip-sample interaction can be calculated from 
the potential energy stored in the capacitor formed be- 
tween the tip and the back-electrode. The (real) force 
acting on the tip F z , which is considered constant dur- 
ing a virtual arbitrary infinitesimal tip displacement 5z, 
performs a virtual work 5w = F z ■ Sz — —SU, where 
U = U c + Ub is the total energy of the system includ- 
ing contributions from both the capacitor and the bias- 
ing battery which maintains a fixed potential difference 
V between the both electrodes. In response to this dis- 
placement, the battery transfers a charge SQ between the 
electrodes in order to keep their potential difference fixed. 



It costs a change of SUb = —SQ ■ V in the energy of the 
battery. Whereas the energy of the capacitor changes by 
SU C — \SQ ■ V, which implies Ub = —2U C , i.e. 

5U = SU C + 6U b = -SU C . 

The electrostatic force is therefore 

_ 5U_ SU C _ ISC 2 
oz oz I oz 

and is always attractive because SC/8z < 0. 



Appendix B: Conducting sphere against a thick 
dielectric slab 



The force between a conducting sphere of radius R 
at potential V facing a dielectric slab grounded on the 
bottom can be calculated by means of the image charge 
method. For a semi-infinite dielectric, we found that the 
solution is given by remarkably simple generalization of 
the treatment in section 5.08 of Smythe's textbook^ 8 - for 
a semi-infinite conductor. Details, further analytic re- 
sults and useful approximations, which are of general in- 
terest for scanning force microscopy, will be presented 
elsewhere 42. The potential 4> in the region between the 
sphere and the slab is generated by a series of point 
charges {q n ,z n } inside the sphere and their correspond- 
ing images {— f$q n , —z n } below the surface of the dielec- 
tric, where = (e — eo)/(e + eo)i e an d £o being the per- 
mittivities of the dielectric and of vacuum, respectively. 
The first charge q\ = 47reoUi? is located at the center 
of the sphere Z\ = R + s. Physically, the image charges 
represents the effect of the polarization induced at the 
surface of the dielectric which causes a jump discontinu- 
ity in E z . Together with the other charges they ensure 
that the sphere surface remains equipotential at V. 

We find 



q n = qisinhaf -— ), 

V smn na I 

where cosh a = Zi/R and 



Z n = -Rsinha cothna, 



(Bl) 



(B2) 



as in Smythe's treatment ((3 — 1). Except at the con- 
tact point (s = 0), the charges q n decay exponentially 
fast towards zero and this solution provides convenient 
expressions for the capacitance C(s) = q sp h/V ', where 



q sp h = J^ q n = 47re -RUsinha x ^ 



00 an-l 



P T ' 



sinh ? 



-, (B3) 



n— 1 n— 1 

and the z-component of the force (dC/ds)V 2 /2 

F M = 27re U 2 Y^ —. ■ (cotha - ncothna),(B4) 

"^ sinn na x ' 

n=2 
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the force gradient dFj^/ds, and the electric field at p 
0,z = 



E z = 



V 1 + /3 ^ /J^sinhna 
R sinh a ^ cosh 2 na 



(B5) 



These series have been used to evaluate the solid lines 
in Figs. [5] and |6fb) The truncation error can be re- 
duced by generalizing a trick proposed for the conducting 
sphere-plane problem^ 4 -. If a series is truncated at some 
n = k, the remainder can be summed up analytically 
if one assumes z n> k ~ z^, = R sinh a. Thus, because 

q n +i/q?i>k ^ /3R/(zi + Zoo), 



n=k+l 



9fe+i 



9fe+i 



1 - £R/(*i + Zoo ) 1 - /3e" 



Adding the correction to the first 5 terms, q sp h is ob- 
tained with an accuracy of 10 -5 for e/e = 5.9, s/R = 
0.1. 

For an ideal conductor ((3 — 1) these expressions re- 
duce to those of Smythe 58 and diverge in the limit s40 
(i.e. a — > 0). In our case /3 < 1, and in the same limit 
the resulting series converge and can in fact be summed 
explicitly ^2- The result for F(s = 0) was given with- 
out proof in Eq. (2) of Ref, 51 . For NaCl (e/e = 5.9, 
/3 = 0.71) one obtains limiting values of C/ireoR = 6.98, 
F/ttcqV 2 = -6.77 (i.e. F = -0.188 nN/V 2 independent 
of sphere radius), F' /ne^R- 1 = 188.7 and E/VR- 1 = 
20.4. In the case of a dielectric slab of finite thick- 
ness t, mirror images of all previously mentioned charges 
with respect to the grounded back-electrode must also 
be considered. They ensure that the field lines inside the 
dielectric become perpendicular to the grounded back- 
electrode instead of spreading radially. If zi = R + s <C t, 
further images charges induced by those mirror charges 
can be neglected to order 0((zi/t) 2 ). Moreover, the elec- 
tric force exerted by the mirror charges on the biased 
sphere can then be approximated as the Coulomb force 
between q spn at its center and a lumped mirror charge 
-(1 - I3)q sph , 2(zi + t) away, i.e. 



Fad 



-(1-13) qi P h T/2 2e 

= Ke V 



A-iver 



(2tf 



R\ 2 /q sp h\ 2 



e + eo V. t J V q 1 J 



For relevant values s < R ~10 nm and t ~ 1 mm, this 
correction is below (R/t) 2 /6.8 ~ 10~ n times the force 



given by Eq. (|B4I) . i.e. negligible in practice. A simi- 
lar expression of comparable magnitude was proposed in 
Ref. 42 , but was erroneously assumed to represent Fm- 



Appendix C: Discretized integrals for finite tip 
oscillation amplitudes 

Assuming that N + 1 equispaced data points {zi\ are 
sufficiently close together such that g(z) remains almost 
constant within an interval length 5 — 2A/N, the inte- 
gration in Eq. (|16p can be approximated by a finite sum 



(g(z)) 



i 



N 

■E 

i=0 



w i9i 



where g % = g(z t ) i s either g % = a(z { ) or g t = C(z % + h): 
Since W(() = 1/^A 2 -C, 2 we obtain 



tv r+ r- 

Wi = W(()dC = arcsin(-^-) - arcsin(-^-) 



*± 



.4 



.4 



where £* = (i ± |)<5 — A are the midpoints between 
Q and Ci±i- Taking into account the rapid variation of 
W(Q near the integration limits defined as (q = —A and 
(^ = A, the square root singularities of W(C) at those 
turning points are approximately included with this mod- 
ified trapezoid integration method. Sufficiently far from 
those points Wj ~ W(Q)5 so that the standard trapezoid 
approximation is recovered. The analogous approxima- 
tion for Eq. flTT)) namely 

1 N 

( g '(z))^-yw: 9t 



-o 



involves^ 



W* 



A 



it 



(W(C)dC 



c 




Note that in the i-)0 limit only the data points at the 
two limits are taken into account. Indeed, if TV = 1, 
A = 8/2 and W = W\, hence (g) = (g + ffjv)/2, 
and Wq = -W*, hence (g 1 ) = (gN - go)/ 2 A, so that 
Eqs. (J14I15I) consistently approximate the correspond- 
ing zero-amplitude equations, Eqs. (|12I13[) . Similarly, 
if N=2, A = 6 and one obtains W = W 2 , W\ = and 
W * = -W£,W? = and Eqs. (J12I13J) are again recov- 
ered. 
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